First-Principles Calculations of Luminescence Spectrum Line Shapes for Defects in 
Semiconductors: The Example of GaN and ZnO 



m 
o 

(N 



Audrius Alkauskas, John L. Lyons, Daniel Steiauf, and Chris G. Van de Walle 

Materials Department, University of California, Santa Barbara, California 93106-5050, USA 

(Dated: March 14, 2013) 

We present a theoretical study of broadening of defect luminescence bands due to vibronic cou- 
pling. Numerical proof is provided for the commonly used assumption that a multi-dimensional 
vibrational problem can be mapped onto an effective one-dimensional configuration coordinate di- 
agram. Our approach is implemented based on density functional theory with a hybrid functional, 
resulting in luminescence lineshapes for important defects in GaN and ZnO that show unprecedented 
agreement with experiment. We find clear trends concerning effective parameters that characterize 
luminescence bands of donor- and acceptor-type defects, thus facilitating their identification. 

PACS numbers: 63.20.kp, 61.72.Bb, 71.55.-i, 78.55.Cr 
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Defects play a key role in the properties of solids. From 
the early days of color centers, the study of luminescence 
and absorption has been crucial to defect characteriza- 
tion [l[ . Theoretical efforts to calculate the broadening of 
optical transitions at defects due to the interactions with 
lattice vibrations were pioneered by Huang and Rhys [2} 
and Pekar J3| . While those theories and their general- 
izations [l|, |j| have been very successful in describing the 
shape of experimental optical bands [l|, l5j , this inevitably 
required the use of empirical fitting parameters. Theory 
has thus been limited in its ability to aid the microscopic 
identification of defects or produce accurate predictions. 

In this Letter we report that unprecedented precision 
can now be achieved by rigorously mapping the multi- 
dimensional vibrational problem onto an effective one- 
dimensional configuration coordinate diagram, combined 
with advanced electronic structure techniques [1, 0] • We 
demonstrate the power of the approach with the example 
of a number of defects in GaN [5| and ZnO [8|] , two tech- 
nologically crucial wide-band-gap semiconductors. Ex- 
cellent agreement with experiment is achieved for well- 
characterized defects, and new insights into vibronic cou- 
pling emerge. 

Our electronic structure calculations are based on den- 
sity functional theory using the hybrid functional of Ref. 
Q in the VASP code The fraction a of the screened 
Fock exchange admixed to the semilocal exchange was set 
to 0.31 for GaN and 0.36 for ZnO to reproduce experi- 
mental band gaps (3.5 eV and 3.4 eV). By describing bulk 
electronic structure better and reducing self-interaction 
errors, hybrid functionals substantially improve the accu- 
racy of defect calculations |10h12| . Defects were treated 
via the supercell approach [6j , the interaction with nuclei 
was described within the projector-augmented wave for- 
malism Q , and electron wave functions were expanded in 
plane waves with a cutoff of 400 eV. Normal modes and 
frequencies have been calculated using finite differences. 

We illustrate the methodology with the example of the 
MgGa acceptor in GaN, a crucial defect since it is the 
only acceptor impurity capable of making the material 
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FIG. 1: (Color online) ID configuration coordinate diagram 
describing optical absorption and emission at a point defect. 
The minima of the ground-state and excited-state potential 
energy surfaces are displaced. AEr e>g ^ are the relaxation en- 
ergies and fi{ e>g } the effective phonon frequencies. ZPL in- 
dicates the zero-phonon line, i.e. the transition between the 
zero-point vibrational states in excited and ground-state con- 
figurations. 



p type. While electrically acting as a shallow impurity 
with modest ionization energy, optically MgQ a behaves 
as a deep center [l3[ : recombination of an electron at the 
conduction-band minimum (CBM) with a hole localized 
on the neutral Mg^ acceptor gives rise to a broad blue 
luminescence band @, [l| . The calculated zero-phonon 
line (ZPL) energy (see Fig. [I} is 3.24 eV. We assume here 
that optical transitions start with a delocalized charge 
carrier; excitonic effects are small [14l |. 

The general theory of luminescence was described in 
Refs. [3, 0, H[. When optical transitions are dipole- 
allowed, as is the case for the defects studied in this 
work, at T=0 the normalized luminescence intensity 
(lineshape) in the leading order (the Franck-Condon ap- 
proximation) can be written as G(hui) = Cu> 3 A(huj), 
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where A(hui) is the normalized spectral function 



Mjfru) = ^2\(x e o\Xgn)\ 2 S(E ZPL -fkv gn -fkj), (1) •§ 



and C _1 = J A(hu)uj 3 d(hoj) . The sum runs over all vi- 
brational levels with energies fvuj gn of the ground state, 
X are ionic wave functions, and Ezpl is the energy of the 
ZPL. Generalization to finite T is straightforward. 

Evaluation of A{TwS) is complicated by the fact that, 
first, the sum includes all relevant vibrational degrees of 
freedom, and second, normal modes Q c and Q g in the 
excited and the ground state are usually not identical. 
The two are related via the Duschinsky transformation 
Q c = JQ g + AQ and (XeolXgn) ar e thus highly- 

multi-dimensional integrals. For small molecular systems 
recursive techniques to calculate such integrals have been 
developed [l6| and implemented [TtJ • The large number 
of vibrational modes that occur for defects in solids ren- 
der such a direct approach computationally prohibitive. 

Broad optical bands have most often been described 
via ID configuration-coordinate diagrams (CCDs) HQ 
[Fig.Q], based on the assumption that the large number of 
vibrational modes (with different frequencies) contribut- 
ing to the lineshape can be replaced by a single effective 
mode (sometimes a small number of modes). The pa- 
rameters entering the ID model are the modal mass M 
of the effective vibration, the displacement of the po- 
tential energy minima AR, and the effective frequencies 
il g and f2 [Fig. Q]. Based on these, the widely used 
"Huang-Rhys (HR) factors" @ are defined as the average 
number of phonons created during a vertical transition: 
S g = AE g /hQ g and 5* c = AE C /Ml g . There are many 
examples where a ID model with empirical fitting pa- 
rameters provides a good approximation to experimental 
luminescence lineshapes 0, Q ; still, because it is strictly 
valid only when all the modes have the same frequency 
0], its general applicability has often been questioned. 
More importantly, the use of fitting parameters precludes 
linking to potentially valuable microscopic information 
about the defect and limits the predictive power. 

Here we address this problem using the following strat- 
egy. Vibrations that couple strongly to the distortion 
of the geometry are expected to be dominant in A(huj). 
Such modes have finite weight on the atoms that experi- 
ence the largest relaxations, i.e., the atoms close to the 
defect. When only a small number of atoms are included, 
(XeolXgn) can be evaluated exactly, taking mode mixing 
into account [1(3, Il7j . This exact evaluation can then 
serve as a test of the accuracy of any approximations. 
Once an approximate treatment has been validated in 
this fashion, it can be applied to much larger systems of 
atoms provided it is sufficiently less numerically demand- 
ing than the exact evaluation. 

In the case of Mgg a a hole is localized on a N neigh- 
bor of the Mg atom, and the five atoms surrounding this 
hole account for more than 90% of the whole relaxation. 
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FIG. 2: (Color online) Normalized luminescence spectrum of 
the MgGa(0/— ) transition, calculated taking into account vi- 
brational modes of only those atoms that relax most (labelled 
in the inset). Black solid line (and shaded area): calculation 
including mode mixing; blue dashed line: parallel-mode ap- 
proximation; red solid line: effective ID vibrational problem. 



The calculated luminescence lineshape, taking mixing be- 
tween the resulting 15 vibrations into account, is shown 
in Fig. [3J Multi-dimensional overlap integrals were cal- 
culated using the molfc code [l7j | . We have applied a 
Gaussian smearing with a small ct=0.01 eV to simulate 
additional broadening mechanisms, resulting in a smooth 
lineshape. 

Recursive algorithms [l6| lead to exploding computa- 
tional requirements when applied to larger atom clusters. 
A simplification that is often used for molecules 17] and 
almost always implied in solids 0] is the parallel-mode 
approximation, in which the eigenmodes of either the 
ground state or the excited state are chosen as common 
vibrational states, leading to J=l. (XeolXgn) then factor- 
izes into ID integrals, each corresponding to one vibra- 
tional mode, greatly reducing computational complexity. 
As seen in Fig. [21 the resulting lineshape is indeed close to 
the exact result. Therefore, while mode mixing is present, 
it is not substantial. 

Now that we have validated the parallel-mode approx- 
imation we can include more atoms, since overlap inte- 
grals become easy to calculate. However, the number of 
terms that have to be included grows very rapidly with 
system size. This reflects the fact that important modes 
often do not occur in the gap of the bulk phonon spec- 
trum but are resonances, and therefore not well localized 
in real space. Consequently, further approximations are 
required. This we achieve by devising a suitable ID CCD 
based on computed parameters as outlined below. 

The weight by which each normal mode k contributes 
to the distortion of the defect geometry during optical 
transition can be written as pk = (AQk/AQ) 2 , where 
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Here a labels atoms, i—{x, y, z}, AR a i = R - a ,i — Rg,a,i 
is the distortion vector, R{ e , g };a,i are atomic coordinates, 
and qk-ai is the unit vector in the direction of the normal 
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mode k (J2ai Qk;aiqi;ai = Sk,i). We find that it is useful 
to define an effective frequency 



O 2 
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^P{c,g};feW { 2 , g};fc , (3) 
fc 



where W{ e ,g};fc is the frequency of the mode 5{ e ,g} ; fe- 

Parameters -Ezpl, AQ, Jl g and Sl c define a ID CCD (cf. 
Fig. [T]) for a quantum oscillator with unit mass and can 
be used to calculate the luminescence lineshape. Gaus- 
sian smearing is still applied, but it should now reflect 
the replacement of many vibrations at various frequen- 
cies with one effective frequency. Inspection of the mean- 
square deviation for the distribution of phonon frequen- 
cies that contribute to the distortion leads to er=0.025 
cV (~ 0.6M7 g ). The result for Mg Ga in Fig.[2]shows that 
our rigorously defined ID model is an excellent approxi- 
mation to the multi-dimensional calculations. 

Now that the validity of the ID model has been estab- 
lished, it can be extended to larger numbers of atoms 
and we no longer need to explicitly calculate normal 
modes and frequencies to determine the effective param- 
eters AQ [Eq. ©] and fy Cig} [Eq. ©]. Indeed, by 
inserting the expression for AQ& into the one for AQ 
in Eq. ([2]), one can show that when all the atoms in 
the supercell are included in the vibrational problem, 
(AQ) 2 = J2 a i TO aAi? 2 [i . The modal mass is defined via 

AQ = M^ 2 AR, where (Ai?) 2 = J2 a ,i AR li- Effective 
frequencies fi can be obtained by mapping the potential 
energy surface around the respective equilibrium geome- 
tries along the path that linearly interpolates between 
the two geometries 18(. A third-order polynomial fit 



was found to suffice in all cases. The frequency SI in the 
quadratic term defined in this way is equivalent to the 
one calculated from Eq. ([3]). Third-order anharmonic 
corrections that affect vibrational wavefuntions and thus 
overlap integrals, were included in the calculations of the 
spectral function perturbatively. For all subsequent cal- 
culations, we have used 96-atom wurtzite supercells, and 
relaxations of all the atoms were included in determin- 
ing effective parameters. The resulting parameters are 
summarized in Table |U 

The luminescence lineshape G(fka) for Mgc a a t T=0 
K is shown in Fig. [31(a), together with ex per imental low- 
temperature data from Refs. [l9[ and (20| . To facil- 
itate comparison with experiment, the calculated line- 
shapes were shifted to bring the maximum of the lumi- 
nescence in agreement with that of the experimentally 
measured curves. The magnitude of the shift provides 
an estimate for the error in the -Ezpl (thermodynamic 
transition level), and is less than 0.1 eV for all defects. 
This error bar reflects both the remaining inaccuracy 
of even the most advanced first-principles methods, and 
any electrostatic corrections (such as excitonic effects or 
donor-acceptor interactions) not included in the present 
model. The width of the theoretical band (0.44 eV) is 
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FIG. 3: (Color online) Calculated (solid lines) and measured 
(symbols) luminescence lineshapes for various defects in GaN 
and ZnO. (a) MgGa in GaN, expt. from Refs. 19l (disks) and 
[U (stars) (b) C N in GaN, expt. from Ref. [H; (c) V N in 
GaN, expt. from Refs [2j|. (d) No in ZnO, expt. from Ref. 
[2q ] The arrows indicate the ZPL, and the calculated spectra 
were shifted by ~0.1 eV, as discussed in the text. 



only slightly larger than the experimental value of 0.36- 
0.37 eV. We derive the effective vibrational frequency in 
the ground state Jl g =47 meV, and the HR factor S s —12. 
The good agreement between the calculated and exper- 
imental lineshapes for MgGa attests to the power of the 
approach presented here, when used in combination with 
state-of-the-art density functional calculations. 

The agreement with experiment is even better for our 
second example, Cn- This defect has been sug gested as 
a source of yellow luminescence (YL) [a, |22j, |23j based on 
the transition C^j- 



•e — >C N [21j (where e is an electron 



at the CBM). This defect again exhibits hole localization 
in the neutral c harg e state, but now the hole is localized 
on the C atom [21(. The calculated luminescence line- 
shape at T=77 K [Fig. [3(b)] agrees very well with the 
measurements of Ref. [23j , in which YL was convincingly 
attributed to a C-related defect. Our calculated effec- 
tive frequencies (M!{ e g }=36, 42 meV) and HR factors 
(5{e, g }=10, 11) are in an excellent agreement with those 
measured in Ref. [H] (Table UJ. 

Our next example is Vn- Nitrogen vacancies have low 
formation energies in p-type GaN, and they have been 
suggested [24| as a cause of the YL in p-GaN observed in 
Ref. [25)]. The calculated luminescence lineshape for the 

[Fig. [3(c)] shows impressive 



transition Kt 3 +e — > 



y+2 



agreement with these low-temperature experiments. 

To demonstrate that the approach is not limited to 
GaN, in Fig. [3(d) we compare the calculated and mea- 
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TABLE I: Effective parameters for various defect-related luminescence transitions in GaN and ZnO. AQ and AR: total mass- 
weighted and total distortions; M: modal mass; fi{ e , g }: effective frequencies in the ground and excited states (charge state 
in parentheses); -Ezpl: zero-phonon line energy; FWHM: full- width at half maximum of the band; T: temperature for which 
the FWHM is given; S{ C:g y. Huang-Rhys factors. If experimental parameters are not explicitly given in the corresponding 
experimental papers, they are extracted using the original data and are shown in italics. 
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sured [26J luminescence lineshapes at T=300 K for the 
deep Nq acceptor in ZnO, corresponding to the transi- 



— >N [L2j. The agreement between theory 



tion N Q +e" 
and experiment is again extremely good. 

Our ability to calculate accurate parameters allows us 
to examine some general trends. For the two substitu- 
tional acceptors in GaN analyzed above, total distortions 
AR amount to 0.22 - 0.24 A, and HR factors to 10 - 12, 
irrespective of whether the hole is bound to C or N. Other 
acceptors with anion-bound holes (Voa, Beca, Zno a ) 
show very similar behavior. The donor Vn, on the other 
hand, is very different. Its defect wave function is com- 
posed mainly of Ga 4s states, and AR associated with 
the (+3/+2) transition is 0.45 A, almost twice as large 
as in the case of acceptors. The distortion mostly affects 
the four nearest Ga atoms, leading to a large modal mass 
and small effective frequencies (Table U). In conjunction 
with large relaxation energies (0.72 and 0.82 eV) this re- 
sults in very large HR factors (<S{ e ,g} = 35, 36). This is in 
contrast with the acceptors, where anions are involved in 
the distortion, leading to a smaller modal masses, higher 
effective vibrational frequencies, and hence smaller HR 
factors. Similar trends are observed for ZnO: we find 
^{e,g} =28-40 meV and S^ c ^ = 15-23 for acceptors with 
anion-localized holes (No, Vz n , Liz n ) while fi.f2{ Cig }=16, 
21 meV and S , { e g }«50 for donors with cation-derived 
states (Vo)- The general result for acceptors in ZnO 
is in accord with experimental data of Ref. [27. 

While such a posteriori interpretations are simple and 
intuitive, they are only reliable if based on an accu- 
rate microscopic description of the defect. We note that 
model calculations have yielded results that were very 
different from our first-principles values (e.g., 5'=3.5— 6.5 
for defects related to YL in Ref. [28j]), starkly illustrating 
the shortcomings of such approaches. 

High values of HR factors mean that it is very difficult 
to determine the ZPL in experimental luminescence spec- 
tra. Indeed, the weight of the ZPL is expotentially su- 



eo|Xgo)| ~ exp{-5 {c!g} } (equal- 
]). As seen in Fig. such com- 



pressed for larger S: \() 
ity holds for 57 g = tt c : 
plication arises for all the defects studied here and is es- 
pecially aparent for Vn- This highlights the practical use 
of calculations exemplified in the current work. 

The examples have demonstrated that our methodol- 
ogy is capable of producing luminescence lineshapes in 
very good agreement with experiment, as well as quan- 
tities that can be directly compared with experimental 
parameters. The achieved agreement is based, of course, 
on the accuracy of the underlying electronic structure 
method, but also on the applicability of the ID model 
to broad luminescence bands. While the explicit consid- 
eration of many vibrational modes is sometimes needed 
to understand the experimental spectra when S w 01 
(Ley defects with moderate electron-phonon coupling) 
[E 0, 0] , We have demonstrated that such a model, with 
suitably calculated parameters, is indeed valid for defects 
with large electron-phonon coupling (S 1), even if 
many phonon modes couple to the optical transition. An- 
other important conclusion that follows from our work is 
that the effective mode frequency is usually much smaller 
than that of LO phonons (91 meV in GaN and 73 meV 
in ZnO), contrary to what has been commonly assumed 
in phenomenological approaches Q . Our conclusion is in 
full agreement with detailed experimental measurements 
for color centers in alkali halides [30] , indicating the gen- 
erality of the obtained result. 

Our findings are important for future studies of semi- 
conductors and insulators that exhibit broad defect lu- 
minescence bands. In particular, as our examples showed 
the methodology provides a means for identifying the mi- 
croscopic origin of the numerous as yet unassigned lumi- 
nescence bands in technologically important wide band- 
gap materials. More generally, our work attests to the 
success of first-principles methods to describe electron- 
phonon interactions in solids f3~H beyond the application 
to perfect crystals. 
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